clear all;
clc;

patnum = 2;
dir1=['patient',num2str(patnum),'_08192014'];
dir2=['p',num2str(patnum),'_act_curves'];
addpath(dir1);
addpath([dir1,'\',dir2]);

% ----------------------- Geometry Parameters -----------------------------

% import data
patientgeom = load('femspine_patientgeom.txt');
cordgeom = load('patient_cordgeom.txt');

epidur  = 0; % 0 = intradural, 1 = epidural
subloc  = 1; % 1 = 1 mm above cord; 2 = 1 mm below dura
th_elec = -10; % 0, -10, -20

% location of electrode

eletag_epi = 'epi_';
fdist_elec_epi = cordgeom(patnum,8);

eletag_sub = 'sub_';
if(subloc==1)
    fdist_elec_sub = cordgeom(patnum,9);
elseif(subloc==2)
    fdist_elec_sub = 1-cordgeom(patnum,9);
else
    return;
end

ptag=['p',num2str(patnum),'_'];
loctag_epi = [ptag,'d',strrep(num2str(fdist_elec_epi),'.','p'),'_'];
loctag_sub = [ptag,'d',strrep(num2str(fdist_elec_sub),'.','p'),'_'];

% ----------------------- Spine Geometry ----------------------------------

% cord geometry
delr_dura = 0.2;

betageom = patientgeom(patnum,:);
a_cord = betageom(1)/2;
b_cord = betageom(2)/2;
a_csf  = betageom(3)/2-delr_dura;
b_csf  = betageom(4)/2-delr_dura;
a_dura = a_csf + delr_dura;
b_dura = b_csf + delr_dura;
ycenter = 13.61;
csf_center = [0,ycenter];
ytopepi = 25;

y_cord = (ycenter-b_csf) + betageom(5)*2*b_csf;
cord_center = [0,y_cord];

% ------------------------ Electrode Location -----------------------------

ytmp = ycenter+b_csf+delr_dura;
yelec_epi = ytmp+(ytopepi-ytmp)*fdist_elec_epi;

yelec_sub = (y_cord+b_cord)+(ycenter+b_csf-(y_cord+b_cord))*fdist_elec_sub;

elec_center_epi = [0,yelec_epi];
elec_center_sub = [0,yelec_sub];

th_rad_0 = 0*(pi/180);
th_rad_n10 = -10*(pi/180);
th_rad_n20 = -20*(pi/180);
Arot_0   = [cos(th_rad_0),-sin(th_rad_0);sin(th_rad_0),cos(th_rad_0)];
Arot_n10 = [cos(th_rad_n10),-sin(th_rad_n10);sin(th_rad_n10),cos(th_rad_n10)];
Arot_n20 = [cos(th_rad_n20),-sin(th_rad_n20);sin(th_rad_n20),cos(th_rad_n20)];

% epidural
rp_epi_0 = Arot_0*[0;yelec_epi-y_cord];
xe_epi_0 = rp_epi_0(1); ye_epi_0 = rp_epi_0(2) + y_cord;
relec_epi_0 = [xe_epi_0,ye_epi_0];

rp_epi_n10 = Arot_n10*[0;yelec_epi-y_cord];
xe_epi_n10 = rp_epi_n10(1); ye_epi_n10 = rp_epi_n10(2) + y_cord;
relec_epi_n10 = [xe_epi_n10,ye_epi_n10];

rp_epi_n20 = Arot_n20*[0;yelec_epi-y_cord];
xe_epi_n20 = rp_epi_n20(1); ye_epi_n20 = rp_epi_n20(2) + y_cord;
relec_epi_n20 = [xe_epi_n20,ye_epi_n20];

% subdural
rp_sub_0 = Arot_0*[0;yelec_sub-y_cord];
xe_sub_0 = rp_sub_0(1); ye_sub_0 = rp_sub_0(2) + y_cord;
relec_sub_0 = [xe_sub_0,ye_sub_0];

rp_sub_n10 = Arot_n10*[0;yelec_sub-y_cord];
xe_sub_n10 = rp_sub_n10(1); ye_sub_n10 = rp_sub_n10(2) + y_cord;
relec_sub_n10 = [xe_sub_n10,ye_sub_n10];

rp_sub_n20 = Arot_n20*[0;yelec_sub-y_cord];
xe_sub_n20 = rp_sub_n20(1); ye_sub_n20 = rp_sub_n20(2) + y_cord;
relec_sub_n20 = [xe_sub_n20,ye_sub_n20];


% white matter
th_unit = 0:pi/20:2*pi;
rcord   = a_cord*b_cord./sqrt((b_cord*cos(th_unit)).^2+(a_cord*sin(th_unit)).^2); 
xcord   = cord_center(1) + rcord.*cos(th_unit);
ycord   = cord_center(2) + rcord.*sin(th_unit);
% csf space
rcsf   = a_csf*b_csf./sqrt((b_csf*cos(th_unit)).^2+(a_csf*sin(th_unit)).^2); 
xcsf   = csf_center(1) + rcsf.*cos(th_unit);
ycsf   = csf_center(2) + rcsf.*sin(th_unit);
% dura
rdura   = a_dura*b_dura./sqrt((b_dura*cos(th_unit)).^2+(a_dura*sin(th_unit)).^2); 
xdura   = csf_center(1) + rdura.*cos(th_unit);
ydura   = csf_center(2) + rdura.*sin(th_unit);

% grey matter
x1 = [0.00 ,0.62 ,1.24 ,1.87 ,2.09 ,2.24 ,1.45 ,2.03 ,2.22 ,2.03 ,1.21 ,0.62 ,0.00];
y1 = [-0.44,-0.12,0.96,2.04,2.06,1.88,0.43,-1.12,-1.62,-2.12,-2.18,-1.61,-0.88];
y1 = y1+cord_center(2);
xgrey = [x1,-fliplr(x1)];
ygrey = [y1,fliplr(y1)];
         
% ----------------------- Importing Data ----------------------------------

p_tag = num2str(patnum);

pol = 1;
noaxons = 200;

if(patnum<3)
    prefix_epi = ['pos_',eletag_epi];
    prefix_sub = ['pos_',eletag_sub];    
else
    prefix_epi = eletag_epi;
    prefix_sub = eletag_sub;
end

if(epidur==1)
    prefix = prefix_epi;
    loctag = loctag_epi;
else
    prefix = prefix_sub;
    loctag = loctag_sub;
end

% DC fibers
dcthresh_LT1 = load([prefix,'LT1p5','_dc_act_',loctag,'t0','.txt']);
[dcthresh_LT1(:,1),dc_si_LT1_0] = sort( dcthresh_LT1(:,1)*pol );

dcthresh_LT2 = load([prefix,'LT6','_dc_act_',loctag,'t0','.txt']);
[dcthresh_LT2(:,1),dc_si_LT2_0] = sort( dcthresh_LT2(:,1)*pol );

dcthresh_TT1 = load([prefix,'TT1','_dc_act_',loctag,'t0','.txt']);
[dcthresh_TT1(:,1),dc_si_TT1_0] = sort( dcthresh_TT1(:,1)*pol );

dcthresh_TT2 = load([prefix,'TT3','_dc_act_',loctag,'t0','.txt']);
[dcthresh_TT2(:,1),dc_si_TT2_0] = sort( dcthresh_TT2(:,1)*pol );

dcthresh_ATP = load([prefix,'atp','_dc_act_',loctag,'t0','.txt']);
[dcthresh_ATP(:,1),dc_si_ATP_0] = sort( dcthresh_ATP(:,1)*pol );

% DR fibers
drthresh_LT1 = load([prefix,'LT1p5','_dr_act_',loctag,'t0','.txt']);
[drthresh_LT1(:,1),dr_si_LT1_0] = sort( drthresh_LT1(:,1)*pol );

drthresh_LT2 = load([prefix,'LT6','_dr_act_',loctag,'t0','.txt']);
[drthresh_LT2(:,1),dr_si_LT2_0] = sort( drthresh_LT2(:,1)*pol );

drthresh_TT1 = load([prefix,'TT1','_dr_act_',loctag,'t0','.txt']);
[drthresh_TT1(:,1),dr_si_TT1_0] = sort( drthresh_TT1(:,1)*pol );

drthresh_TT2 = load([prefix,'TT3','_dr_act_',loctag,'t0','.txt']);
[drthresh_TT2(:,1),dr_si_TT2_0] = sort( drthresh_TT2(:,1)*pol );

drthresh_ATP = load([prefix,'atp','_dr_act_',loctag,'t0','.txt']);
[drthresh_ATP(:,1),dr_si_ATP_0] = sort( drthresh_ATP(:,1)*pol );

% DC location
dc_xy = load([ptag,'dc_xyzloc.txt']);
dc_xy = dc_xy(:,1:2);
dc_xy = dc_xy + ( cord_center'*ones(1,noaxons) )';

dc_xy_LT1_0 = dc_xy(dc_si_LT1_0,:);
dc_xy_LT2_0 = dc_xy(dc_si_LT2_0,:);
dc_xy_TT1_0 = dc_xy(dc_si_TT1_0,:);
dc_xy_TT2_0 = dc_xy(dc_si_TT2_0,:);
dc_xy_ATP_0 = dc_xy(dc_si_ATP_0,:);

% ----------------------- DC Activated before DR --------------------------

nocase = 5;
xy_selact = cell(1,nocase);
relec = zeros(nocase,2);

nodc_LT1_0 = sum( dcthresh_LT1(:,1) < drthresh_LT1(1,1) );
nodc_LT2_0 = sum( dcthresh_LT2(:,1) < drthresh_LT2(1,1) );
nodc_TT1_0 = sum( dcthresh_TT1(:,1) < drthresh_TT1(1,1) );
nodc_TT2_0 = sum( dcthresh_TT2(:,1) < drthresh_TT2(1,1) );
nodc_ATP_0 = sum( dcthresh_ATP(:,1) < drthresh_ATP(1,1) );

xy_selact{1} = dc_xy_LT1_0(1:nodc_LT1_0,:);
xy_selact{2} = dc_xy_LT2_0(1:nodc_LT2_0,:);
xy_selact{3} = dc_xy_TT1_0(1:nodc_TT1_0,:);
xy_selact{4} = dc_xy_TT2_0(1:nodc_TT2_0,:);
xy_selact{5} = dc_xy_ATP_0(1:nodc_ATP_0,:);

rp_sub_0 = (cord_center + rp_sub_0')';
relec = (rp_sub_0*ones(1,nocase))';

title_graph = cell(1,nocase);
title_graph{1} = 'LT-1.5';
title_graph{2} = 'LT-6';
title_graph{3} = 'TT-1';
title_graph{4} = 'TT-3';
title_graph{5} = 'AT';

% ----------------------- Spatial Activation ------------------------------

figure;
for k = 1:nocase
    
    perdc = 100*length(xy_selact{k})/noaxons;

    subplot(2,3,k);
    hold on;
    fill(xcord,ycord,[1,1,1],'LineWidth',3);
    plot((xcsf+xdura)/2,(ycsf+ydura)/2,'k--','LineWidth',4);
    fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
    plot(relec(k,1),relec(k,2),'ko','MarkerSize',10,...
         'MarkerFaceColor','w','LineWidth',4);
    plot(xy_selact{k}(:,1),xy_selact{k}(:,2),'k.','MarkerSize',25);
    set(gca,'XTick',[],'YTick',[]);
    
    % p1(y)=8.5, p2(y)=18.5
    text(-3,18.5,[num2str(perdc,2),' %'],'FontSize',36);
    title(title_graph{k},'FontSize',36);
    axis tight;
    axis square;
    hold off;
    
end

figure;
subplot(2,3,k);
hold on;
fill(xcord,ycord,[1,1,1],'LineWidth',3);
plot((xcsf+xdura)/2,(ycsf+ydura)/2,'k--','LineWidth',4);
fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
plot(relec(k,1),relec(k,2),'ko','MarkerSize',12,...
     'MarkerFaceColor','w','LineWidth',4);
plot(xy_selact{k}(:,1),xy_selact{k}(:,2),'k.','MarkerSize',25);
set(gca,'XTick',[],'YTick',[]);
axis off;

L1 = 'white matter'; 
L2 = 'dura mater'; 
L3 = 'grey matter'; 
L4 = 'electrode location';
L5 = 'DC fibers activated';
legend(L1,L2,L3,L4,L5,'Location','BO');
set(gca,'FontSize',36);